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We suggest that Hamiltonian matrix elements between physical states in QCD might be approximated, in 
Coulomb gauge, by “lattice-improved” tree diagrams; i.e. tree diagram contributions with dressed ghost, trans¬ 
verse gluon, and Coulomb propagators obtained from lattice simulations. Such matrix elements can be applied 
to a truncated-basis treatment of hadronic states which include constituent gluons. As an illustration, we apply 
this hybrid approach to the heavy quark potential, for quark-antiquark separations up to 2.4 fm. The Coulomb 
string tension in SU(3) gauge theory is about a factor of four times greater than the asymptotic string tension. 
In our approach we show that a single constituent gluon is in principle sufficient, up to 2.4 fm, to reduce this 
overshoot by the factor required. The static potential remains linear, although the precise value of the string 
tension depends on details of the Couloumb gauge ghost and gluon propagators in the infrared regime. In this 
connection we present new lattice results for the transverse gluon propagator in position space. 


I. INTRODUCTION 

It is not obvious that particle-like gluons, which are vi¬ 
tal to perturbative QCD, really make sense as constituents of 
hadrons, particularly highly excited hadrons with higher spin. 
Perhaps such states can only be described by strings of some 
kind which connect to quarks. The picture of linear Regge tra¬ 
jectories as arising from a spinning line-like object is certainly 
compelling, and, if there is any sense in which gluons are con¬ 
stituents of hadrons, then surely the hrst challenge is to hnd 
out whether this line-like object has a substructure that can be 
understood in terms of individual gluons. The simplest case 
to study is the lowest energy state containing a static quark- 
antiquark pair separated by a distance R. The line-like object 
connecting the quarks must manifest itself as a color electric 
flux tube. Does this flux tube have a substructure that involves 
individual gluons, as in the “gluon chain” proposal of ref. [1]? 

To address this question, and in fact to even define what is 
meant by an individual gluon, it is necessary to work in a fixed 
gauge. We will use Coulomb gauge, which has the advan¬ 
tage that a confining potential is already built into the dressed 
Coulomb propagator. This fact has been verified numerically 
in many lattice studies [2-7] and it is, moreover, a necessary 
condition for a non-vanishing asymptotic string tension [8]. 
The problem, however, is that the SU(3) Coulomb string ten¬ 
sion Ocoui derived from the instantaneous Coulomb propaga¬ 
tor is a factor of four times larger than the asymptotic string 
tension [2, 3], which seems too much of a good thing. We 
may ask whether constituent gluons in a static quark-antiquark 
state can somehow reduce the Coulomb string tension to the 
known asymptotic value. If that turns out to be true then 
we could go on to study other hadronic states, such as the 
low-lying glueballs, which would be another natural setting in 
which to investigate gluons as constituent particles. 

Among the physical states in Coulomb gauge, containing 


e.g. a static quark and antiquark pair, are superpositions of 
states of the form 

f " 

•' i=\ 

q\^)Aki{xi)Ak^{x2) ■ ■ ■Ak,XXn)q\R)\Q)vs.c , ( 1 ) 

where |0)vac is the true vacuum state, q\q^ are massive quark- 
antiquark creation operators, and 'P is a function which, in 
a variational approach, may depend on some set of param¬ 
eters. Suppose we have a finite, not necessarily orthogonal, 
set of such states, labeled by an integer (Ij)}. From these a 
set {|j)} of orthogonal states can be constructed. If we could 
compute Hamiltonian matrix elements {j\H\k) in the Hilbert 
space spanned by this truncated basis, then the standard proce¬ 
dure is to diagonalize the Hamiltonian in the truncated basis, 
minimize the energy of the lowest energy state by adjusting 
the variational parameters, and in this way arrive at an esti¬ 
mate for the static quark potential at separation R. A similar 
strategy could be employed in spectrum calculations, involv¬ 
ing states with dynamical quarks. The problem, of course, 
is to calculate the relevant overlaps and Hamiltonian matrix 
elements. In principle this task can be carried out by lat¬ 
tice Monte Carlo in Coulomb gauge, and this was the path 
followed in [9]. The problem was also addressed from the 
Dyson-Schwinger point of view in [10]. In this article we will 
suggest a somewhat different approach, inspired by renormal¬ 
ized perturbation theory. 

Consider the expression 

Cjk{t) = {j\e-^Xk) ■ ( 2 ) 

The time derivative evaluated at f = 0 gives us {i\H\k), while 
Cjk{Q) is the overlap {j\k). The prescription for calculating 
Cjk{t), in ordinary perturbation theory, is essentially the same 
as the prescription for calculating an S-matrix element: Sum 


2 


all of the tree diagrams which contribute to this expression, in¬ 
cluding all n-point vertices. The vertices are the one-particle 
irreducible (IPI) n-point functions appearing in the quantum 
effective action, and the propagators are the full (or “dressed”) 
propagators of the theory. The task of perturbation theory is 
to compute these renormalized propagators and IPl n-point 
functions. This program can be carried out analytically, within 
the limits of validity of an asymptotic expansion, if the spa¬ 
tial separations of the quarks and gluons in states Hy)} are 
small. If this is not the case then the program fails, because the 
perturbative expansion for the relevant propagators and ver¬ 
tices rapidly diverges. Let us observe, however, that only a 
finite number of tree diagrams contribute to the calculation of 
Cjk{t). The sum of trees is not an infinite expansion; it is fi¬ 
nite and, given the dressed propagators and n-point vertices, 
it supplies the exact answer, regardless of the magnitude of 
the renormalized coupling. Therefore, if propagators and rel¬ 
evant vertices could be calculated by some non-perturbative 
approach, say by Dyson-Schwinger equations or lattice Monte 
Carlo simulations, then the tree diagrams could be summed, 
and Cjt[t) could be calculated. From those quantities, the 
spectrum of the Hamiltonian in the subspace of Hilbert space 
spanned by the set {\j)} could be calculated. 

In this article we will take a first step along these lines, 
by taking ghost, transverse gluon, and Coulomb propagators 
from lattice Monte Carlo simulations, neglecting all vertices 
apart from those arising from the non-polynomial operator in 
the Coulomb gauge Hamiltonian. We will use the resulting 
tree diagrams to compute Hamiltonian matrix elements in a 
truncated basis of thirteen states, consisting of one state with 
no constituent gluons, and twelve states with a single con¬ 
stituent gluon in various spatial distributions. It is found that 
the static quark potential derived in this way remains linear, 
but the asymptotic string tension depends on an overall con¬ 
stant factor associated with the ghost propagator. In the ab¬ 
sence of decisive data on this point, we simply tune the factor 
in the ghost propagator to get the known result. Hopefully 
future lattice Monte Carlo investigations of the ghost propa¬ 
gator will make our study more predictive, but for now we 
only show that inclusion of constituent gluons in static quark- 
antiquark states, along the lines of the gluon chain model, pro¬ 
vides a very plausible mechanism for reducing the string ten¬ 
sion from the pure Coulombic value, which is much too high, 
to the value consistent with numerical simulations. 

Below in Section II we review some of the numerical re¬ 
sults and conjectures in ref. [2], which motivate the work pre¬ 
sented here. In Section III we will present our proposal in 
detail, and in particular explain how the non-local operator 
which appears in the Coulomb gauge Hamiltonian is treated 
in the tree diagram approach. Expressions for the Hamilto¬ 
nian matrix elements in terms of lattice-improved tree dia¬ 
grams are derived in section IV, and they require, in addi¬ 
tion to the Coulomb propagator already obtained in [2], also 
the ghost and transverse gluon propagators. In Section V we 
will show our lattice Monte Carlo results for the equal times 
transverse gluon propagator in position space, relying on refs. 
[5, 11, 12] for the infrared behavior of the ghost propagator. 
In Section VI we will bring all these results together, com¬ 


pute an estimate for the static quark potential, and show how 
the superposition of zero and one constituent gluon states can 
bring the static quark potential down from the Coulomb result 
by a large numerical factor dependent on the ghost propaga¬ 
tor. We conclude in Section VII. Finally, in the Appendix we 
discuss the results of Monte Carlo simulations of gluon, ghost 
and Coulomb propagators in the context of Dyson-Schwinger 
approach. 


II. COULOMB POTENTIAL AND THE GLUON CHAIN 
MODEL 

In a recent article [2] we calculated the non-perturbative 
Coulomb potential in SU(3) pure gauge theory via lattice 
Monte Carlo simulations in Coulomb gauge, and found it to 
be 

Vcoui{R) = (ycouiR “ ’ ( 3 ) 

where the Coulomb string tension is 

Ocoui « 20.5(4) fm-2 

Ri (893±9MeV)2 , (4) 

(see also [3]) which is about four times the accepted value of 
the asymptotic string tension a = (440 MeV)^. This is clearly 
too much of a good thing. While it is helpful that a linearly 
confining potential is obtained in Coulomb gauge by what 
amounts to one (dressed) gluon exchange, it is not acceptable 
that the static quark potential is four times too large. Moreover 
there is no indication that the color electric field generated by 
the Coulomb propagator around the quark-antiquark sources 
is restricted to a flux tube. Something else important must be 
going on. 

The key point is that Vcoui{R) is not necessarily the minimal 
energy of a state containing a pair of static quark-antiquark 
sources. It is, rather, the interaction energy of a specific state 
in Coulomb gauge, namely 

|0)w=^^(0)^^(/?)|0)vac , (5) 

where 

fl>o[A] = (A|0)vac (6) 

is the true vacuum wave functional. The energy of this 
state, including self-energies, is given by the logarithmic time 
derivative of the Euclidean-time correlator 

V{R) = -lim^log5,(0|c-^'||0)5, . (7) 

As Dirac indices and quark kinetic energies are not relevant 
to our study, it is sufficient to compute, in a Euclidean action 
formulation, the logarithmic time derivative of a correlator of 
short timelike Wilson lines 

V{R) = -lim|log(Tr[L,(0)Lj(R)]) , (8) 
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where 


Lt (jt) = T exp 


ig [ dtAo{x,t) 

Jo 


(9) 


In a lattice formulation, the calculation of V(R) boils down to 
to calculating the logarithm of the vacuum expectation value 
(VEV) of products of two timelike links, evaluated at equal 
times. This method for computing the Coulomb energy was 
first suggested in [7], and the potential defined in this way 
includes an /^-independent self-energy term. The self-energy 
term, proportional to the inverse lattice spacing, can be iden¬ 
tified and subtracted away, with the result for the //-dependent 
quark-antiquark interaction potential Vcoui{R) shown above in 
(3) and (4). The coefficient 7r/12 in (3) is the correct value for 
the Liischer term, but this coefficient is not simply assumed. 
Rather, it appears to be the likely continuum limit of the val¬ 
ues derived at finite lattice spacings, cf. [2]. 

Since the state lO)^^ cannot be the minimal energy state 
containing static quarks, it is reasonable to consider states 
with n constituent gluons of the form shown in (1), where 
what we mean by n “constituent gluons” is simply that there 
are n A-field operators that operate on the true ground state. 
In diagramatic representation this ket vector corresponds to n 
transverse gluon lines emerging from a blob at time fo = 0, and 
the bra vector is a blob with n transverse gluons lines entering 
at some later time t. The static quark and antiquark lines at¬ 
tach to either end of the blob. The gluon chain model of ref. 
[ 1 ] proposed that the minimal energy state containing a static 
quark antiquark pair consists of some number n of constituent 
gluons arranged roughly in a line (or cylindrical region) be¬ 
tween the quark and antiquark, with the color ordering of the 
gluons correlated with the spatial ordering of the gluons along 
the line. The original motivation was the observation that the 
force between colored sources, in two-loop perturbation the¬ 
ory, grows very rapidly with R when the running coupling ap¬ 
proaches values of 0(1). So the idea is that as a quark and 
antiquark separate, the energy grows until at some point it is 
energetically favorable to place a gluon in between the quark 
and antiquark, which halves the effective color charge separa¬ 
tion. As the quark and antiquark continue to separate, even¬ 
tually it is energetically favorable to introduce a second gluon 
between the two sources, and so on. In the end, the mini¬ 
mal energy state would contain approximately n—R/Ro con¬ 
stituent gluons, where Rq is essentially the average distance 
between gluons, and the kinetic plus inter-gluon interaction 
energy is Eq. Then the total energy of the long chain is approx¬ 
imately E{R) = hEq = {Eo/Ro)R, which is a linear potential 
with string tension Eq/Rq. 

It now seems clear that this picture is untenable, because it 
assumes that at some point, as color charges separate, the po¬ 
tential between colored sources grows faster than linear. But 
if the Coulomb potential between gluons is instead asymptot¬ 
ically linear, then there is no advantage to introducing more 
gluons. If the gluons were arranged exactly in a line between 
the quarks, and the inter-gluon Coulomb potential is linear 
with string tension a', then the overall Coulomb energy would 
be g'R regardless of the number of gluons in the chain. The 


kinetic energy of the gluons, and the transverse fluctuations 
away from the line joining the quark-antiquark pair, could 
only increase this energy. So it would appear that the mini¬ 
mal energy really is the zero constituent gluon state (5), and 
we have found that this state has a string tension which is far 
above the asymptotic string tension. 

However, this conclusion ignores the fact that a state with 
a fixed number of constituent gluons is not an eigenstate of 
the Hamiltonian, and there will always be non-zero matrix 
elements {m\H\n) between states with different numbers of 
constituent gluons. This means that the lowest energy state 
is certain to be a superposition of states with different num¬ 
bers of constituent gluons. Diagonal matrix elements contain 
both purely kinetic contributions, indicated schematically in 
Fig. 1(a), and instantaneous Coulomb interactions, shown in 
Fig. 1(b). Matrix elements between states with different num¬ 
bers of constituent gluons are associated with diagrams such 
as Fig. 1(c). All matrix elements are derived from the time 
derivatives of the diagrams indicated. In ref. [2] we showed 
that if we make some plausible assumptions about the behav¬ 
ior of the diagonal and off-diagonal Hamiltonian matrix ele¬ 
ments, then the static quark potential associated with the low¬ 
est energy state could be reduced from the purely Coulomb 
value by a large numerical factor, while retaining the asymp¬ 
totic linearity of the potential. For details of the model cal¬ 
culation we refer the reader to [2]. In the present article we 
will actually evaluate the type of diagrams just indicated, in 
a “lattice-improved” tree-diagram framework, to see if they 
really do have the conjectured effect. 


III. LATTICE IMPROVED TREE DIAGRAMS EOR T HE 
STATIC QUARK POTENTIAL 

A. Preliminaries 

For completeness and to establish notation we begin with 
the usual preliminaries regarding Coulomb gauge. The 
Coulomb gauge Hamiltonian is // = Hgi^e + Hcoui + Hmatter, 
where 

Hgiue = -B"), 

Hcoui = y d\d^y /-^-p%x)/K-\x,y-A)p‘’{y)/-K 

( 10 ) 

with 

p‘‘ = p^+pI+pY 

Ji = -SI.9{A) , ^=det[.^]. (11) 

Here p|(.r) = gq]{x)tfjqj{x), p|(x) = gqi{x)tfjq]{x) and 
Pg(x) = —gf‘‘^‘^Ag{x)Ej^{x) are the charge density of quarks, 
antiquarks and gluons, respectively, and ^^(A) is the covari¬ 
ant derivative. H„,alter is the part of the Hamiltonian contain¬ 
ing dynamical matter fields. It will not be needed here, since 
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FIG. 1. Diagrams which, after a time derivative, contribute to Hamiltonian matrix elements, (a) The graph which determines the kinetic energy 
of constituent gluons, (b) One of the graphs determining the Coulomb energy of an u-gluon state. The blob labeled “K” is the instantaneous 
non-polynomial Coulomb operator, (c) Schematic of a graph contributing to an off-diagonal Hamiltonian matrix element between states with 
n and n + 1 gluons. Here one of the A-field operators in the Coulomb operator K{A) contracts with a gluon in the final state. 


in this article we are only concerned with static color sources, 
which can be represented by Wilson lines in the time direc¬ 
tion. The operator-ordering terms ^*2 do not appear at the 
classical level, and therefore do not appear in the construction 
of tree diagrams. 


and the /T-propagator 



= . (16) 


An equivalent Euclidean path-integral formulation in first- 
order formalism [13-16] can be based on the generating func¬ 
tional 


Z[J] = 

j DA*l J DEY exp j d^x {lEiAi - ^ {e} + B ■) - /. 
j dtd'^xd^y{pg + gJ4)xjK[x,y,t,A]{pg + J4)y^t 


— iJA, 


( 12 ) 


The Coulomb propagator demands some special attention; 
it is not the same as the fT-propagator. It is defined as 

Dt(x-y) 


= 

\z 5J^{x)5J^(y) J 
= {K‘“’{x-y-A))S{x4-y4) 

+ J dhdMK“^ix-z)p^g{z)p^g{w)K^'\w-y)) 


= 5‘>‘’{K{x-y)5{x4-y4)+Pix-y) 


(17) 


where the “f r” superscript in the measure indicates that the ind = 4 dimensions P(x - y) may have both an instantaneous 
integration is restricted to transverse A,£ fields, and Pg(x) = non-instantaneous part 

Pix-y)=Pi„^(x-y)5{x4-y4)+P„on{x-y) , (18) 


We define the transverse gluon propagator 

D‘ff{x,t) = {Anx,t)A‘’m) 

= Dij{x,t)5‘“’, (13) 

the E —A propagator 

&;f{x,t) = {Ef{x,t)A']m) 

= 5,y(x,f)5“\ (14) 

the ghost propagator 



= G(x-y)5“\ (15) 


so that 

D‘Yl{x-y) = 8“^(^K(x-y)+Pinst{x-y)^5{x4-y4) 

+Pnon{x-y). (19) 

The instantaneous part of the Coulomb propagator will be de¬ 
noted 

K‘‘\x-y) = 5‘‘‘’K{x-y) 

= 8‘^‘’(^K{x-y)+Pi„st{x-y)^ . (20) 

It was shown by Zwanziger [13, 15] that both g^D'Y^{x — y) 
and g^K{x — y) are renormalization group invariants. For now 
only the instantaneous part of the Coulomb propagator will 
be needed. We will denote the relationship of K to operator 
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expectation values by 


static quark 


K={K) + {KppK)i„„ , (21) 


where the subscript on the last term indicates the instanta¬ 
neous part of the VEV. The relationship between the Coulomb 
potential between static sources in the fundamental represen¬ 
tation and the instantaneous Coulomb propagator is 

Vcoui{R) = -g^CFK{R) . ( 22 ) 

For later use we will introduce the notation 

Kc{R) =Vcout{R) 

Gc{R)=gG{R). (23) 


Our graphical notation for these propagators is shown in 
Fig. 2. We will assume that the transverse gluon propagator 
has the form 


Di 


, , f / kikj 

~ J {27tr ki+coi ■ 'w 


j 


(2n:)4 kl + col 
d^k 

(2;r)3 2(0k ' 


k^ 

k 


kikj 


= 5ijD{x,t)-D\f{x,t) 
1 

2 ' 


D{x,t) = -5ijDij{x,t) , 


(where k^ = k- k), and also that 


Di 


f d^k 

k^e'^tiot+kx) 

(3.. 

J {27t:Y 

kl + (ol 



= 5ijD{x,t)-D‘i]\x,t) 


D{x,t) = ^5ijDij{x,t) . 


(24) 


(25) 


At the perturbative level (0^ = |fc|. At the non-perturbative 
level (Ok will be something else, with 


(Ok = 


massive propagator 
y/W + jk"^ Gribov propagator 


(26) 


as possible candidates. Our approach is to determine the rel¬ 
evant gluon propagators in position space from lattice Monte 
Carlo simulations. We will not need to know cOk explicitly, 
although we will, in Section V, compare our numerical results 
for the transverse gluon propagator in position space with the 
forms implied by the above candidates for (Ok- However, the 
assumption that Dij^Dij have the form above allows us to con¬ 
clude that 

D' (x, 0) = lim (x,t) 
t=o dt 

= -\8\x) 

limD(x,e) = ^sign(e)5^(x) . (27) 

e->0 2 


_ Coulomb 

. ghost 

^VWWVW transverse gluon 

propagator 

FIG. 2. Diagrammatic notation for Coulomb gauge propagators. 


B. Decomposition of the Coulomb vertex 


In the tree diagram framework we suggest here, there is a 
question of how the non-polynomial operator 

p%x)K‘‘\x-y,A)p‘’{y) (28) 

should be handled. Fet us first consider the case where the 
charge operators p contract only with “external” operators in 
the initial and final states, either heavy quarks or constituent 
gluons. These operators will be denoted p'^”. Then we also 
have to consider the possibility that some A-field operators 
in the perturbative expansion of K‘‘^{x — y,A) contract with 
gluon operators in the initial and final states. In some dia¬ 
grams there are no such contractions, and those diagrams sum 
up to the dressed /T-propagator (K). We must also consider 
the product of pKp operators in which two charge operators 
contract with initial and final states, i.e. p‘^^KppKp‘^^. Again 
there are diagrams in which there are no contractions of oper¬ 
ators in either K with the external constituent gluons. Those 
diagrams sum up to (KppK), and add to the /T-propagator 
(K) to produce the dressed Coulomb propagator/r(x — y) (see 
( 21 )). 

In other diagrams, however, it is necessary to consider a 
“Coulomb vertex,” as indicated schematically in Fig. 1(c), 
where one or more field operators in K contract with operators 
in the initial and/or final states. It turns out, as we will now 
show, that Coulomb vertices can be decomposed into products 
of ghost operators. Coulomb operators, and transverse gluon 
operators that contract with the external states. In our tree 
diagram formulation, these ghost operators and Coulomb op¬ 
erators just become dressed ghost and Coulomb propagators, 
and this prescription amounts to a partial resummation of the 
full perturbation series. 

We will begin with the case so it is sufficient to 

just consider the expansion of the Coulomb operator 

<(A) = G-(-V2)G^^, (29) 

where the ghost operator is 

{-V-DY^ = , 


(30) 
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The perturbative expansion of the ghost operator begins with The Coulomb operator is then 


- (JC - y) + (31) 


or 


(32) 

This equation can be solved iteratively, and from here on we 
will drop both color and spatial indices. The solution is a 
power series 


where we define 


= -4;^ E 

( ^ ) m=0n=0 


(35) 


In the perturbative expansion of K, each of the operators 
M contains a single A operator, and in a Coulomb vertex we 
have to choose one or more of these to connect to external 
gluons in the initial or final states. Let us begin with the case 
of a Coulomb vertex with a single gluon emerging. Denote by 
M* the operator which contains the A-field contracting with an 
external gluon field, and the resulting operator (corresponding 
to the blob with one gluon line coming out) we’ll denote K'^. 
Then we have 


M = gfA,di 


1 


1 OO N-l 

(34) E + 1) E ' 

^ A'^=0 m=0 

1 °° “ 

= 7z^E E("*+«+2)™*^" 

I / m=0/7=0 

. OO OO 1 

= E^E4'« + '^+2)^'”(^/A*5)p^M" 

1 OO 1 

= E^(“+E"" 

+p^E""(s/'’'‘>)p^E(»+i)«" 

= K{gfA*d)G + G{gfA*d)K. (36) 


Now let’s consider a blob with two gluon lines coming out. 
By the same reasoning 


1 OO N-2 N-m-2 

= E (^+ E E 

V ) N=0 m=0 n=0 

= E E f^{m + n + k + 3)M"‘M*M"M*M'^ 

V ^ ) m=0n=0k=0 

= p^E^EE(”+"+*+3)""W-4‘'^)p^""(«/''‘'^)p^"‘ 

= E^<"+E«"feM‘^)p^ E"‘ 

T OO 1 OO ^ OO 

+p^E__«"'(s/^’<>)p^E("+i)«"(s/-<*<^)p^E«* 

-I OO 1 OO -I OO 

= K{gfA*d)G{gfA*d)G + G{gfA*d)K{gfA*d)G + G{gfA*d)G{gfA*d)K . 


( 37 ) 
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The general case follows from induction, and is indicated 
schematically in Fig. 3. The tree diagram approximation is 
to replace the operators G by the corresponding propaga¬ 
tors, as shown in the figure. In diagrammatic terms this means 
we either neglect completely contractions of field operators in 
one {K or G) blob with field operators in another blob, or we 
assume that the only effect of such contractions is to multiply 
the point-like vertex by a constant, which can be absorbed into 
a multiplicative factor c in the ghost propagator. 

At this stage it would appear that the Coulomb vertices in¬ 
volve /T-propagators rather than Coulomb propagators. How¬ 
ever, this neglects the fact that there is a second contribution to 
vertices with N external gluon lines, which comes about from 
expanding the KppK operator. We recall that the Coulomb 
propagator is actually a sum of the K-propagator and {KppK). 
A similar statement is true for Coulomb vertices. Let us con¬ 
sider one term in the expansion of the K operator with N 
external gluon operators, in which there are n products of 
G{gfA*d) to the left of the K operator, and N — n products 
of {gfA*d)G operators to the right of the K operator, i.e. 

i — II—y i — N — n—y 

Ti = G{gfA*d )... G{gfA*d) K {gfA*d)G... {gfA*d)G . 

(38) 

Again, in the tree diagrammatic expansion we neglect any di¬ 
agrams which connect G operators to K operators or to other 
G operators, and the result is a product of ghost propagators, 
A'-propagators, and A field operators which contract with op¬ 
erators in the initial and/or final states. To each term of this 
kind, there is a corresponding term which arises from the ex¬ 
pansion of the KppK operator: 

72 = I G{gfA*dY^G{gfA*d) TTpj 

r ^N-n^ '] 

X |p7: igfA*d)G... {gfA*d)G j . (39) 

In the tree diagram expansion the G operators become ghost 
propagators, and KppK becomes (KppK). Then in the tree- 
level approximation to Ti + T 2 the K operators appear in the 
combination {K) -f {KppK), which is simply the Coulomb 
propagator K. This is completely general. The effect of 
adding the KppK operator in computing Coulomb vertices 
in this framework is equivalent to considering only the expan¬ 
sion of the K operator, and then replacing the TT-propagator 
with the Coulomb propagator. 

One might ask why we have neglected terms in the expan¬ 
sion of KppK such as 

I^G{gfA*d)... G{gfA*d)K{gfA*d)G ... {gfA*d)Gp | 

X ^PG{gfA*d) ■ ■ ■ G{gfA*d)K{gfA*d)G... (g/A*5)G| . 

(40) 

The answer is that all of the diagrams which contribute to the 
VEV of such a quantity involve loops which cannot be ab¬ 
sorbed into a ghost or Coulomb propagator. In the tree dia¬ 





p 





p 


(tree diagrams) cj) 



P 



P 



FIG. 3. Decomposition of a vertex, generated by the Coulomb op¬ 
erator K{x — y,A), into products of operators. In the tree diagram 
decomposition, G and K operators become ghost and Coulomb prop¬ 
agators respectively, as explained in the text. 


gram framework we propose here, contributions of that kind 
are dropped. 


IV. THE HAMILTONIAN IN A BASIS OF ZERO AND ONE 
CONSTITUENT GLUONS 


In the toy model of a gluon chain, explained in ref. [2], the 
average number of gluons in the chain will grow linearly with 
the quark-antiquark separation. If this idea really works, then 
there should then be some intermediate distance range, pre¬ 
sumably just after the onset of confinement, where the average 
number of gluons is less than one, i.e. the minimal energy state 
is mostly a superposition of just zero and one-gluon states. We 
will denote a finite set of such states, containing also a static 
quark-antiquark pair of separation/?, as {|n)}, where |0) is the 
zero-gluon state 


|0)=^t“(0)^^“(/?)|0)v.., (41) 


and the remaining one-gluon states (u > 1) have the form 


«)=?'’'“ (0) 


J d^x^•\"\x)Af{x)t^l^ 


q^"{R)\0)v 


(42) 


with ta = the SU(3) group generators. From these states 
we can construct an orthonormal basis, and diagonalize the 
Hamiltonian in that basis. For this purpose we need to com¬ 
pute 


Hmn = {m\H\n) , Omn = {m\n) 


(43) 
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The trial one-gluon wavefunctions may depend on 

some variational parameters, which are chosen such that the 
lowest energy state obtained by diagonalization is minimized. 

There is a big simplification if we choose our trial one- 
gluon wavefunctions (x) to be transverse, V • = 0. 

Then, since the transverse gluon propagator will always con¬ 
tract with the index of an external wavefunction, and the co¬ 
ordinate X of the external wavefunction is always integrated 
over, we can simply drop the kikj/k^ piece of the transverse 
projection operator, because that piece will always act like a 
divergence on Effectively, then, we are allowed to drop 
D[j^(x,f) in (24), and just use 


Dij{x,t) dijD{x,t) , 
d^k 1 


Dixd) = j 

Likewise, for e « 0, 


{iTty 2(Ok 


e<kxe-(0kt 


(44) 


Dij{x,e) ^5;;sign(e)5^(x) . (45) 


{k) defined by 

(^) = / ^ ■ (47) 

Then 

(«|n) = 3Cf j d^x f*^"\x)fl"\x) . (48) 

So (x) is the analog of an ordinary quantum mechanics 
one-particle wavefunction, and ideally we would like to first 
choose fi"\x) based on some combination of intuition and 
convenience, transform to ff'\k), and then compute 

in) 

from (47). But in practice it would be difficult to derive 'Pj 
analytically, because of the ^/Wk factor in (47). Even numeri¬ 
cally this would be challenging, because the integrand of (47) 
is an oscillating function. It is simpler to work directly with 

A. The truncated basis 


We also denote, at equal-times, D{x) = D{x,t =0). With 
transverse 'Pj ^ (x) and « > 1 

Omn = {m\n) 

= Txtata / d^Z\d^Z2'i‘f'^*{z\)D{z\-Z2)'i‘f\z2) ■ 

(46) 

with Tr{tata) = 3 Cf where Cf is the quadratic Casimir in the 
fundamental representation. 

Now, for any choice of transverse 'p|"^(x), there exists 
another transverse function fj'"\x), with Eourier transform 


Eor the numerical work in this article, we will choose as an 
ansatz transverse wavefunctions of the form 


»p(") = V X 



Fap{x,y,z) 


-xdzFap 

-ydzFap 

2Fap +xdjcFap +ydyFap 


with a = 1,2,..,A, p = 1,2,..,M and n = {a 
Then Fy^p is chosen to have the form 


(49) 
l)M-f j3. 


Fap{x,y,z) = fa{z)L'p_l 



(50) 


where r = \J 2 ? is the polar coordinate in the transverse 

direction, is an associated Laguerre polynomial, a is a vari¬ 
ational parameter, and 

1 a = 1 

. (51) 

The motivation for this choice of F^p is the fact that if a gluon 
is located on the z-axis between the quarks, with x = y = 0 
and 0 < z < R, then the linear piece of the Coulomb potential 
between the gluon, the quark, and the antiquark, sums up to 
(JcZ + (7c{R — z) = <7cR which is independent of z. So if the 
gluon is not too far away from the z-axis, it can move more 
or less freely in the z-direction between the quarks. Hence the 
excitations in the z-direction might be as in a one-dimensional 


square well of width Ri R. Since there is no reason for the 
wavefunction to drop to zero exactly at z = Q,R, we allow for 
a “well” a little bigger than R, namely R + 2zo, where in this 
work we have taken zo = a/3. Note that for r » \z\,\R-z\ 

F{x,y,z) ^ exp(-2r/a) . (52) 

The set of orthogonal functions in the integration measure 
yg-^r/a ^j.g j-jjg associated Laguerre polynomials L]„{A-r/a). 

It should be understood that an orthogonal set of F^^p does 
not produce an orthogonal set of one gluon states {]«)}; it 
only specifies a set of linearly independent states which to¬ 
gether span a subspace of Hilbert space. The object is to find 
the lowest eigenvalue of the Hamiltonian in this subspace. If 
we define a matrix [FI] with elements H^n and a corresponding 
matrix [O] of overlaps Omn in the non-orthogonal set of states 
{|n)}, as in (43), the eigenvalues of the Hamiltonian operator 
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restricted to this subspace are obtained by solving the gener¬ 
alized eigenvalue equation 

[//]mW (53) 

This is essentially equivalent to the alternate procedure of con¬ 
structing a set of orthonormal states {|v,j)} from the {|n)} 
via, e.g., the Gram-Schmidt procedure, computing matrix el¬ 
ements H'-j = {vi\H\vj) in the orthogonal basis from the H„,„, 
and then diagonalizing the matrix H' in the usual way. 


tained by displacing the initial and hnal state in time, so that 
the transverse gluon propagator is at unequal times, and then 
taking minus the time derivative of the diagram shown in Fig. 

5: 

=-lim^TrfaG j d^zid^Z2^f”\zi)D{zi-Z2,t)^'\"\z2) 

= -TTtJaj dhid\2 'i‘f^\zi)D'{Zl -Z2,0)'f|"'(Z2) 

= 3Cf Jdh]^^f\z)¥t\z) . (55) 

From (47), this could also be expressed in terms of ft as 


B. Matrix Elements 

We now list the improved-tree contributions to each of the 
required Hamiltonian matrix elements H„,„ of eq. (43). The 
state with a static quark-antiquark pair and zero constituent 
gluons will be denoted |0), and states with the static quarks 
and one consituent gluon will be denoted |n) with n > 0. The 
expression for overlap matrix elements (9„,„ of one-gluon 
states was given in (46), while for the zero-gluon state 
O()o = 3 (a sum over three colors), and Oon = 0. 


1. {0\H\0) 


The matrix element is given by the diagram in Fig. 4, and 
this is the simplest contribution. It is 

Hoo = -g^TT{tata)K{R) = 3KciR) . (54) 


9L 


-9'a 


Id^kc0kf*{k)f,ik) , (56) 

although this is actually not the most useful form if the vari¬ 
ational wavefunction is given in terms of 'F, rather than /). 



FIG. 5. Diagram associated with constituent gluon kinetic energy. 


3- (»l^|0) 


The relevant diagrams are shown in Fig. 6, and we find 


H„Q = -2i^ j d^xd^z'i‘\"\z)D{x-z) 

X {Gc(x)diKc{R - x) + Kc{x)diGc{R - Jc)} 

(57) 


FIG. 4. Diagram responsible for the instantaneous Coulomb poten¬ 
tial. 

2. {m\H\n) kinetic term 

The next simplest term is associated with the kinetic energy 
of the one-gluon state for the case m = n, and that is best ob- 

I 


4. {m\H\n) Coulomb term 


The planar diagrams are shown in Fig. 7. There are two dia¬ 
grams in which the Coulomb propagator ends on the left hand 
side quark, and likewise two diagrams where the Coulomb 
propagator ends on the right hand side antiquark. The sum of 
all contributions is 


HZ' = ^jd\xd \2 'f|”’>(Zi)Z)(zi -Z 2 )¥;'\z 2 ){Kc{Zl)+Kc{R-Zl)+Kc{Z 2 )+Kc{R-Z 2 )] . (58) 


There is also a non-planar contribution to the Coulomb en- be 
ergy, which is shown in Fig. 8. This contribution comes out to 


HZ = ~\Kc{R)0„n 


(59) 
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FIG. 6. Diagram responsible for the Hamiltonian matrix ele¬ 
ment. 



9ta 




-9ta 


FIG. 7. Planar diagrams contributing to the instantaneous Coulomb 
energy of the one constituent gluon state. 



FIG. 8. Non-planar contribution to the instantaneous Coulomb en¬ 
ergy of the one constituent gluon state. 


5. {m\H\n) KGG terms 


The diagrams are shown in Figs. 9(a) and 9(b). “KGG per¬ 
mutations” means permutations in the order of the instanta¬ 
neous Coulomb (K) and ghost (G) propagators. The possi¬ 
ble orders are KGG, GKG, and GGK. With the diagrammatic 
rules as before (trace of matrix generators is clockwise), it 
turns out that diagram 9(b) is simply the complex conjugate 
of diagram 9(a). The resulting contribution is 


ttKGG _ 
^mn 


Cf. 

x\kc {x) dfGc (x - y) Gc{y-R)+Gc {x)dfKc {x - y) d]Gc {y-R)+Gc {x)dfGc (x-y)d]Kc{y-R) 


j cPxcP'ycPucP’v D(u)D{v) | {x — (j — v) + {x — {y — v) 


(60) 


V. THE EQUAL-TIMES TRANSVERSE GLUON The Coulomb propagator follows from (22) and the Coulomb 

PROPAGATOR potential determined previously in [2]. In this section we 

present our lattice Monte Carlo results for D(x —y); the in- 

In order to actually compute the //„„ matrix elements we 
need the full propagators Kc{x — y), Gc{x — y), and D{x — y). 
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+ KGG permutations 
(a) 



+ KGG permutations 

(b) 

FIG. 9. Diagrams contributing to the matrix element. These 

entail a 12-dimensional integration. 


frared behavior of the ghost propagator will be taken from 
refs. [5, 11, 12]. 


A. The question of gauge copies 

Since our calculation is in Coulomb gauge, we will begin 
with some remarks about Gribov copies. It is well known that 
any point in field space where the quantity 

^=^LLR^Tr[t/,(x,f)], (61) 

evaluated on each timeslice, is stationary (V 3 is the number of 
sites on the timeslice), will satisfy the Coulomb gauge condi¬ 
tion 


3 

^{Uk{x,t) — Uk{x—k,t)—h.c. = 0 (62) 

k=l 

at all sites x, which is just the lattice version of V • A = 0. 
Gribov copies are the configurations on a gauge orbit which 
satisfy this condition, and those copies which are local max¬ 
ima of R are said to be inside the Gribov region. The global 
maximum of R is said to be in the “fundamental modular re¬ 
gion.” In practice it is impossible to find the global maxi¬ 
mum of the gauge fixing condition, although some studies use 
simulated annealing in an attempt to find a better gauge copy 
than what might be obtained from, e.g., the ordinary over¬ 
relaxation technique. Here we should remark that in a funda¬ 
mental sense there is no such thing as a “better” gauge copy. 
It is obvious that if we are calculating a gauge-invariant quan¬ 


tity, then the choice of gauge copy is irrelevant. Of course, 
it may be that gauge dependent quantities such as propaga¬ 
tors and vertices will vary somewhat from copy to copy, and 
the expectation value for such quantities could very well de¬ 
pend on exactly how the gauge copy is chosen. On the other 
hand, if one could put together propagators and vertices so 
as to compute a gauge-invariant quantity (such as a scatter¬ 
ing amplitude in ordinary perturbation theory, or the energy 
of some physical state), then it is expected that the gauge vari¬ 
ance of the propagators and vertices should cancel out, and the 
details of Coulomb gauge fixing should be irrelevant. For this 
to happen, however, it is important that the same version of 
Coulomb gauge fixing is used for all gauge-dependent quan¬ 
tities. 

Any lattice Monte Carlo calculation in Coulomb or Landau 
gauge uses some version of what might be called “computer 
gauge.” Lattice configurations are generated according to the 
gauge-invariant probability weighting, and then some deter¬ 
ministic procedure is applied to arrive at a particular Gribov 
copy inside the Gribov region. The procedure we have used 
in [ 2 ], and that we will continue to use here, is the method 
of Fourier acceleration, introduced in [17]. Although we will 
not use simulated annealing to try to obtain “better” gauge 
copies, it may still be of interest to compare values of R in 
Coulomb gauge obtained from simulated annealing followed 
by Fourier acceleration, with the values obtained from Fourier 
acceleration alone. We have carried out these procedures on 
the f = 1 timeslice of 25 lattices of volume 24^, generated at 
gauge coupling j3 = 6 . 0 , with each lattice separated by 2000 
update sweeps. Each timeslice was first fixed to Coulomb 
gauge by the Fourier acceleration method (which stops af¬ 
ter a certain convergence criterion is met), and the value of R 
(which is approximately 0.88 at this value of ]3) was recorded. 
The timeslice was then subject to a random gauge transforma¬ 
tion followed by 1000 simulated annealing sweeps, with the 
final fixing carried out via Fourier acceleration, and the pro¬ 
cedure again stops when the convergence criterion is met. In 
every case, the values of R obtained for each lattice by the 
two procedures differed at most at the fourth non-zero digit, 
i.e. a difference on the order of lO^"*^. A variety of cooling 
schedules and final temperatures were tried out for the simu¬ 
lated annealing steps, and seemed to make little difference to 
this result. It may be that Fourier acceleration is already very 
efficient at maximizing R. 


B. Renormalization 

Propagators computed in lattice simulations are bare propa¬ 
gators. The relation of the bare to the renormalized transverse 
gluon propagator is given by the usual relation 

Z)^7-(x)=z2A7(x), (63) 

where is the wavefunction renormalization factor, which 
depends on both the cutoff and the renormalization scheme. 
Relations between various renormalization constants in 
Coulomb gauge were worked out in [13], where we find that 
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g^K{x) in terms of bare coupling and propagator equals the 
same expression in terms of the renormalized coupling and 
propagator. From this reference one can also deduce that the 
combination gG{x)Ai^ is renormalization invariant, in which 
case 


Sbare^bareiA) —* ( 6 '^) 

It is worth noting that Coulomb vertices, in Hamiltonian ma¬ 
trix elements, are always multiplied by a factor of g^. This 
means that the tree level Coulomb vertices discussed in sub¬ 
section involve only renormalization group invariant combi¬ 
nations g^K(x — y) and gG(x — y)A(y). 

This observation has an important consequence; it means 
that when states are normalized (i.e. by division by \/Onn), 
the Hamiltonian matrix elements we compute are renormal¬ 
ization group invariant. Therefore if we compute ghost, gluon, 
and Coulomb propagators at a particular lattice coupling )3, 
and use these to calculate the energy spectrum of the static 
quark-antiquark pair plus constituent gluons, the continuum 
limit obtained by extrapolation to j5 will be the same as if 
we had carried through the spectrum calculation using instead 
the corresponding renormalized propagators and vertices. 


C. Numerical results 


We define the space components of A-field on the lattice, in 
Coulomb gauge, as 


Ai{x,t) 


1 

2 iga 


{Ui{x,t)-Ul{x,t)) , 


(65) 


where g = i/b/jS and a = fl(j3) is the lattice spacing in fermis 
at coupling j3, obtained from the Necco-Sommer formula [18] 


expectation value 


^{R) = ^-5u{TT[Ai{x,t)Ajiy,t)]) , (67) 


where R ~ |jt — y |. *. In terms of the propagator 


Df^{R) = {Afix,t)A‘>jiy,t)]) 

= 5‘‘^Dij{R) ( 68 ) 

we have 

^{R) = lduy"‘’D‘lfiR) = ^duDuiR) 

= Id{R) . (69) 


The simulations have been carried out for lattice volumes 
12^16^20^24^28^30^ and j3 = 5.7,5.8,5.9,6.0; there is 
roughly a factor of two between lattice spacings at the low¬ 
est and highest values of p in this range. Figure 10 is a set 
of logarithmic plots which displays the dependence of D{R) 
on lattice volume, at each value of j3. What we see in these 
figures is that the data at the different volumes may have con¬ 
verged at R < 0.6 fm or so, but at larger separations this is 
not the case. At j3 = 5.7,5.8 the data has dropped below the 
x-axis and gone negative in some region. At P ~ 5.9,6.0 the 
data points remain positive, but there is reason to suspect that 
at still higher volumes this data might also go negative in the 
long distance regime. The data is shown on a linear scale in 
Fig. 11. 

It is interesting to see whether the data in the intermediate 
regime, where we may be seeing scaling, can be fit by a stan¬ 
dard ansatz. Define the equal-time propagators corresponding 
to a massive or Gribov form 


and 



d^k e'*-' 
( 2 ;r)^ 2 {k'^ + 


m 

4-7t^R 


Ki (mR) 


(70) 


a{P) = (0.5 fm) exp [- 1.6804- 1.7331 (j3-6) 

-f0.7849(j3 - 6)^-0.4428(j3-6)^] . (66) 



d^k 

(27r)^ 2{k^ +m‘^/k^)^ 


(71) 


We then use lattice Monte Carlo simulation (with Coulomb The latter integral can be expressed in terms of MeijerG func- 

gauge fixing as described above) to compute the equal-times tions; 


^ Note that a mid-link prescription for position would make no difference to 
separation R, given the i — j restriction. 
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gluon propagator, |3=5.9 



gluon propagator, |3=5.8 



gluon propagator, p=6.0 



FIG. 10. Lattice Monte Carlo results, in physical units, for the equal-times transverse gluon propagator at lattice volumes 12^ — 30^ and 
couplings 5.7 < jS < 6.0. Data is shown on a logarithmic scale. 
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(72) 


Since the data has not convincingly converged at R > 0.6 
fm, we can at least try to fit the data in the range 0 < R < 
0.6 fm, where we do seem to see both scaling and volume 
convergence, to either the massive or Gribov forms, i.e. 


D{R) 


cD'^ass{R,m) 

cD(^rih{R^m) 


(73) 


where c and m are the fitting parameters. Fig. 12(a) shows a 
best fit for the massive propagator, and Fig. 12(b) for the Gri¬ 
bov propagator, on a logarithmic scale. The data points shown 
are combined data for all couplings on the largest 30^ lattice 
volume. Clearly neither of these fits are very convincing. The 


fitting functions can get the short distance behavior correctly, 
where the propagator is large, but then they deviate in the 
longer range tail, where the propagator is small. By fitting 
instead to the logarithm of the data one can get quite a good 
fit to the tail, but then the fits both go wrong at the short dis¬ 
tance end, especially when plotted on a linear scale. Neither 
the massive nor the Gribov propagator gives a good account 
of the data in the full 0 <R< 0.6 fm interval. This motivates a 
search for some other functional form for the propagator. We 
have found that this form: 

( 74 ) 
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gluon propagator, |3=5.7 


gluon propagator, p=5.8 
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FIG. 11. Same as Fig. 10, but on a linear scale. 

Gluon Propagator, mas,sive fit Gluon Propagator, Gribov fit 




R R 

(a) massive (b) Gribov 


FIG. 12. Comparison of lattice data for the transverse gluon propagator to fits by (a) a propagator with a mass pole; (b) the Gribov propagator. 


gives an excellent fit to the data. In the logarithmic plot in 
Fig. 13(a), and the linear plot in Fig. 13(b) we show the data 
for D{R) at all four (5 values at the largest 30^ lattice volume. 
On the same plots we also display the best fit (solid line) to 
the data in the interval 0.09 < R <0.6 fm by D‘’'“‘‘{R). The 


constants which give this fit are 

a = 2.35(23) fm-^ , Zi = 5.65(31) fm^^ , c = 0.0469(12) . 

(75) 

There is no particular theoretical justification for the form 
(74). The ansatz just happens to work very well in this particu¬ 
lar distance interval, where we have convergent data. Presum¬ 
ably (74) is simply a good approximation in this interval to the 
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gluon propagator, 30"* lattice 



gluon propagator, 30'* lattice 



R (fm) 


(b) lineal' scale 


FIG. 13. Comparison of the best fit (74) to the equal-times transverse 
propagator data, at various lattice couplings, on the largest 30^ lattice 
volume, (a) logarithmic scale; (b) linear scale. 


bare gluon propagator for this range of couplings is fit by 


D{R) = 0.0469 


exp[—5.65/?^ + 2.35/?] 
^2 


(76) 


where R is in fm and D{R) in fm 2. 

It should noted that there have been previous studies of 
both the transverse propagator in Coulomb gauge, carried out 
in both SU(2) [5, 11, 19] and SU(3) [20] pure gauge the¬ 
ory. These studies computed gluon propagators in momen¬ 
tum space, rather than position space, and the scaling analysis 
was on time asymmetric lattices. For these and other technical 
reasons it is not straightforward to compare our position space 
results directly with the earlier momentum space studies, and 
we will not attempt this here. However, given the strong sen¬ 
sitivity of our results to lattice volume, we believe it would 
be very helpful to carry out further studies of the equal-times 
transverse gluon propagator in both position and momentum 
space, on much larger lattices than those used here. 


D. The ghost propagator 


We now consider the (bare, unrenormalized) ghost propa¬ 
gator 

= (77) 

where the Faddeev-Popov operator on any given timeslice of 
the lattice is 


true, and no doubt very complicated, transverse equal-times 
gluon propagator, which probably violates positivity at larger 
distances, and displays logarithmic corrections at shorter dis¬ 
tances. 


3 r 


^ax,by ' 


:ReTr^ {f,t‘’}{Uk{x) + Uk{x-k))d, 


'xy 


k=]. 

u.bu.a 


-2t‘’fUk{x)5^+^y - ItVUkix - 


(78) 


It is important to note the striking fact, seen in Figures 13(a) 
and 13(b), that the data at different values of the lattice spac¬ 
ing, when expressed in physical units, completely overlap up 
to 0.6 fm or so, beyond which the data still seems to be vol¬ 
ume dependent. This is not expected. In lattice simulations 
we are dealing with the unrenormalized propagator, so data 
sets at different values of the cutoff ought to differ from one 
another by multiplicative constants, and in the coupling range 
shown the largest and smallest lattice spacings differ by al¬ 
most a factor of two. Yet there seems to be no trace of a 
cutoff-dependent multiplicative constant in the range of lat¬ 
tice couplings we have probed. The conclusion appears to be 
that within this range of lattice cutoffs we may have Ri 1, 
at least in the force renormalization scheme which underlies 
the Necco-Sommer formula (66), and if that is the case then 
the bare and renormalized propagators agree in this coupling 
range. Of course must differ very much from unity at suf¬ 
ficiently large j3 and small lattice spacings. In any event, the 


with the 14 (jc) being the spacelike link variables on that times¬ 
lice. However, ^ has eight eigenvectors with zero eigenval¬ 
ues (c = 1 — 8) 

Wax=-^^ac with ^ax.byYbf (79) 

associated with a remnant global SU(3) color symmetry, and 
is therefore not invertible as it stands. We must therefore in¬ 
vert ^ on a subspace orthogonal to these zero modes. An 
alternate procedure is to invert the Faddeev-Popov operator in 
momentum space, as in [5, 11, 12]. Nakagawa et al. [12] find 
in momentum space, for the SU(3) group 

(infrared) (80) 

Burgio et al. [5, 11] and Langfeld and Moyaerts [19] find 
nearly the same power dependence in SU(2). Translating this 
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into position space, we write 


G,(R)=gG(R) = y^^ (81) 

The motivation for multiplying the ghost propagator by cou¬ 
pling g is that, as we have noted above, the combination 
gG{R)Z^^, where gG is the product of the bare lattice cou¬ 
pling and ghost propagator, is equal to the product of the 
renormalized coupling and ghost propagator. We have al¬ 
ready seen evidence that Z^i Ri 1 in the range of couplings 
P € 5.7 —6.0. We therefore treat c(j3) Ricas a constant in this 
coupling range, to be fixed by fitting to some physical quan¬ 
tity. As mentioned below (37), this constant may also absorb 
a contribution correcting for the factorization of the operator 
expectation values shown in Fig. 3. 


VI. THE STATIC QUARK POTENTIAL 

We now define the matrix [H] to have elements H,nn = 
{m\H\n) where Hqq and //„o = were defined in eqs. (54) 
and (57) respectively, and 

IT _ _L ijcoid I Tjnp I ttKGG 

^nm ^nm ' ^nm ' ^nm' ^nm ’ 

where expressions for , 7/,““^, Hnm , have been given 

in eqs. (55), (58), (59), (60). All of these matrix elements, as 
well as the elements C>„,„ of the overlap matrix [O] defined 
in (46) depend on the parameter a that appears in the one- 
gluon wavefunctions ‘F, as well as the parameter c in the ghost 
propagator. 

The procedure, for any given ghost parameter c and quark- 
antiquark separation R, is to compute the energy spectrum of 
the Hamiltonian in a truncated basis which spans the Hilbert 
subspace containing the non-orthogonal zero and one-gluon 
states {!«)}, and vary the parameter a to minimize the lowest 
energy eigenvalue. As we have already mentioned, a possi¬ 
ble procedure is to construct an orthonormal basis from the 
{|u)} states via Gram-Schmidt orthogonalization, compute 
the Hamiltonian matrix elements in this basis, and then di¬ 
agonalize the Hamiltonian. An equivalent method is to solve 
the generalized Hamiltonian equation this 

gives the same spectrum of energy eigenvalues. This proce¬ 
dure is carried out for a range of quark-antiquark separations 
R, and parameters c. If the lowest energy Emin {R) is linear as a 
function of R, then we choose c to get a = Emin {R)/R equal to 
the asymptotic string tension. In our calculation we have cho¬ 
sen indices a,j3 in the wavefunction (see (50)) in the range 
a = 1,2,3,4, P = 1,2,3, so counting the zero-gluon state this 
is thirteen states in all. 

Our final result, and the main result of this paper, is shown 
in Fig. 14. Here we plot the variational/truncated basis es¬ 
timate for the static quark potential V{R) up to /? = 2.4 fm, 
in intervals of 0.2 fm, as compared to the Coulomb potential 
of eq. (3). We have found that as we increase the number of 
states in the truncated basis from two to thirteen, the ground 
state energy converges at around six or seven states. Error 
bars, comparable to symbol size, are estimates of the error de¬ 
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FIG. 14. The Coulomb and gluon chain potentials. 


riving from numerical integrations. The solid line through the 
computed points is a best fit by the function 

V^^R) = af'‘R + 8. (83) 


and we find 


(j/* =4.91(9)fm^^ 

= (437(4) MeV)2 . (84) 

which can be compared to the accepted value of (7 = 
(440 MeV)^. Of course, the parameter in the ghost propagator 
c = 3.5 was chosen to obtain this result, but what is non-trivial 
is the fact that the potential remains linear, despite deviating 
very strongly from the non-perturbative Coulomb potential, 
which we recall had a string tension four times greater than 
the asymptotic string tension. 


VII. CONCLUSIONS 


We have shown how the inclusion of a single constituent 
gluon in the static quark-antiquark state can preserve the linear 
dependence of the energy on quark separation R, at least up to 
R = 2 A fm, while drastically reducing the string tension from 
the pure Coulomb value of (Jcoui = (891 MeV)^. To show that 
the reduction is by precisely the right amount, namely a factor 
of four, will require further numerical studies of the ghost and 
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transverse gluon propagators. 

The next step in our program is to apply our variational 
approach to gluelumps, glueballs, and heavy (but not static) 
quark-antiquark bound states. For quarkonium states the dia¬ 
grams are essentially the same as those shown in section IV, 
only replacing the static quark lines with dynamical quarks 
and antiquarks, and of course there will be wavefunctions for 
relative positions and spin of quark, antiquark, and (zero or 
one) constituent gluons. Likewise, for the glueball states, the 
static quark lines are replaced by gluon propagators. For these 
states, where no constituent is located at a fixed position, it is 
preferable to carry out the calculation of Hamiltonian matrix 
elements in momentum rather than position space. This takes 
advantage of overall momentum conservation, and greatly 
reduces the dimensionality of integrals that must be carried 
out numerically. We hope to report on work along these lines 
in a future publication. 
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the ghost dressing function behaves as 

£/(/?) oc l/p(i+/^)/2 ^ (A.4) 
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For any j3 and ^7(/r) < dc, d{Q) is finite. However, there is 
a range of low momenta where d{p) follows the power-law 
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Appendix: Comparison with Dyson-Schwinger equations 


We note that if D{R) > 0 at all distances, ft)(0) is finite, while 
an IR enhanced co with j3 > 0 implies that the gluon propa¬ 
gator is positivity violating. The momentum space Coulomb 
potential is proportional to the square of the ghost propagator 
and the Coulomb form factor /, 


Cpg^Kip) 


d^{p)f{p) 

n 2 


which satisfies the following DS equation. 


(A.7) 


In this appendix we compare our results for the ghost and 
Coulomb propagators derived from lattice simulations with 
the predictions of Dyson-Schwinger (DS) equations [21-24]. 
In the one-loop, rainbow-ladder approximation, the renormal¬ 
ized DS equation for the momentum space ghost dressing 
function d{p) = gG{p)p^ is given by 


f{p)=m+if[f{p)]-if[m]. 


with 


IfW)] =Ncl 


d^k \ — [p-k)^ d^{k)f(k) 
(2;r)^ 2 o}{k — p) (p — k)^ 


(A.8) 


(A.9) 


1 

d{p) 


^^-I,[d{k)]+Id[d{^)], 


It follows from (A.8) that f{p) is determined up to an overall 
normalization. 


where 


ld[dik)] =Ncj 


d^k \ — (p-k)^ d{k) 
{2n)^ 2co{k — p) (p — k)^^ 


(A.2) 


and 2 (o{p) — D^^{p) is the inverse momentum space gluon 
propagator. Given the function (0{p), the solution of (A.l) 
depends on the value of d(p) at an arbitrary renormalization 
point, which we choose as /i = 1 fm^^ Specifically, the be¬ 
havior of the ghost dressing function in the IR, i.e. for p ^0, 
depends on the value of d{p) and the IR behavior of Co{p). If 
©(O) is finite then so is d{Q), providing that d{p) is less then 
some critical value dc which depends on (O. As d{p) increases 
towards the critical value the ghost dressing function becomes 
IR enhanced, and d{Q) becomes infinite for d{p) = dc with 
d{p)o.\lp^l^ as p —0. In general, if Co(p) behaves in the IR 


A good approximation to the momentum space lattice prop¬ 
agator can be obtained using a simple form motivated by the 
Gribov ansatz 

(0{p) = [p/a)^ + m^(mg/p)^P . (A. 10) 

The functional form adopted in (76), i.e. a positive D(R), 
implies (5=0. In this case the DS equation (A.l) results 
in a ghost propagator which at large R behaves as gG{R) = 
(y/6/8)AR^*/^ withA Ri 1.5. The power behavior is close to 
the lattice fit (81), but the magnitude is approximately a factor 
of two smaller than what is required (according to the analysis 
in this article) to fit the asymptotic string tension. It follows 
from (A.l) that to increase the value of the ghost propagator 
it is necessary to enhance (O in the IR, and this in turn implies 
positivity violation. The lattice data on the gluon propagator 
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does indeed indicate positivity violation. While the analytical 
form in (76) is acceptable in the analysis of the gluon chain 


spectrum, comparison with the DS equation requires taking 
positivity violation into account. 
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